# load data
source("HNC_metadata_tidy.R")
source("data_handling.R")
source("risk_factor_regressions.R")

ancestry <- read.csv("Supplementary_Table_14_data_ancestry.csv") %>% select(-country)

CNburden <- read.csv("Input_mutational_burden/Input_CN_burden.csv") %>% rename(CNV_burden = CN.burden)

cnCOSMIC <- read.csv("Input_signature_attributions/Input_CN_COSMIC_attributions.csv", stringsAsFactors = T, row.names = 1) %>% tibble::rownames_to_column("donor_id")

additional_df <- list(CNburden, cnCOSMIC)

cosmic <- additional_df[[1]]
for (n in 2:length(additional_df)) {
  cosmic <- cosmic %>% left_join(additional_df[[n]], by = "donor_id")
}

# merge with cosmic sigantures and tidy epi variables
data <- data %>%
  left_join(cosmic, by = "donor_id") %>%
  left_join(ancestry, by = "donor_id") %>%
  mutate(
    subsite = factor(subsite, levels = c("Larynx", "OC", "OPC", "Hypopharynx")),
    age_group = factor(age_group, levels = c("55-65", "0-45", "45-55", "65-75", "75+")),
    age = as.numeric(scale(age_diag))
  )

confounders <- c("age", "subsite", "sex", "tobacco_ever", "alcohol_ever")
var <- "African"

result <- NULL
for (v in var) {
  for (sigtype in grep("CN", names(data), value = T)) {
    # Linear regression
    dat <- data[!is.na(data[, var]), ]
    myformula <- paste0(sigtype, " ~ ", paste(c(v, confounders), collapse = " + "))
    print(myformula)
    model <- lm(myformula, data = dat)
    r <- summary(model)$coefficients
    r <- r %>%
      as.data.frame() %>%
      tibble::rownames_to_column("independent_vars") %>%
      select(-`t value`) %>%
      filter(grepl(v, independent_vars)) %>%
      rename(pval = `Pr(>|t|)`) %>%
      mutate(p_adj = ifelse(pval * length(signatures) < 1, pval * length(signatures), 1), .after = "pval") %>%
      mutate(signature = sigtype, .after = "independent_vars")
    result <- rbind(result, r)
  }
}

write.csv(result, "output/Supplementary_Note_Table_7.csv", row.names = F)
